UNIFORMLY HYPERBOLIC ATTRACTOR OF THE 
SMALE- WILLIAMS TYPE FOR A POINCARE MAP IN THE 
KUZNETSOV SYSTEM 



DANIEL WILCZAK 

To the memory of my mathematics teacher, 
Maria Rog 

Abstract. We propose a general algorithm for computer assisted verification 
of uniform hyperbolicity for maps which exhibit a robust attractor. 

The method has been successfully applied to a Poincare map for a system 
of coupled non-autonomous van der Pol oscillators. The model equation has 
been proposed by Kuznetsov T<] and the attractor seems to be of the Smale- 
Williams type. 



1. Introduction. 

Hyperbolic systems of dissipative type, contracting the phase space volume, ma- 
nifest robust attractors. That means that the dynamics does not change qualita- 
tively when the parameters of the system vary in some range. 

There are several models of maps which produce hyperbolic nontrivial attractor 
- like the Plykin attractor or the Smale solenoid [KH) . But for long time there 
were no examples of continuous systems that apparently have hyperbolic strange 
attractors. 

Recently, Kuznetsov [K] proposed a continuous model that exhibits an uniformly 
hyperbolic attractor of the Smale- Williams type in the Poincare map. This is a non- 
autonomous system in and it is constructed on a basis of two coupled van der 
Pol oscillators: 

X = LOqU, 

ii = — wpa; + (Acos(27rt/r) — x^) w + (e/wo)j/cos(ti;ot), 

y = 2a;ou, 

V ^ -2ujQy + (-Acos(27rt/r) ~ y'^) v + {e/2LUo)x'^. 

Let P denote the Poincare map of the above system defined naturally as the shift 
along the trajectories over the period of the vector field 

(2) P{x, u, y, v) = (x(T), u(r), y{T),viT)) . 

Kuznetsov and Sataev |KSj gave a deep numerical study of this system and observed 
that for some range of parameter values there is an absorbing domain for P of the 
toroidal shape as presented in Figure [T] left panel. This domain is a product of a 
3D ball and a circle. 
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Figure 1. Left: an absorbing domain projected onto three coor- 
dinates after a linear change of coordinates. Right: local unstable 
and stable foliations. It is a numerical indication that they are not 
tangent over the absorbing domain. These figures are [KSl Figures 
1 and 7] - used here with the permission of prof. Kuznetsov. 



The observed attractor is of the Smale- Williams type. Numerical studies [Kl 
IKS[ IKSe] gave a strong numerical evidence of the existence of uniformly hyperbolic 
attractor for some range of parameter values around 

(3) wo = 27r, A = 5, T = 6, £ = 0.5 

- see Figure [T] right panel, where the unstable and stable foliations are presented. 

In this paper we propose a method for computer assisted verification that a map 
possesses an uniformly hyperbolic attractor. As a test case we apply the method 
to the model map proposed by Smale [KHj . This is a map defined on the two- 
dimensional solid torus and it is analytically proved to be uniformly hyperbolic. 

Although the method is general, the main motivation for us to undertake this 
study was to prove Kuznetsov's conjecture about the hyperbolicity of the system 
([T]). The following theorem is the main result of this paper. 

Theorem 1. Consider the system with the parameter values Let P de- 
note the Poincare map for this system as defined in Q). There exists a compact, 
connected and explicitly given set B such that 

(1) B is positive invariant with respect to P, i.e. P{B) C B, 

(2) P is uniformly hyperbolic on the maximal invariant set A — f]^^QP^{B) 
with one positive and three negative Lyapunov exponents. 

After a linear change of coordinates the set S is a union of 7 970 392 four- 
dimensional explicitly given cubes. The set yB is a very narrow enclosure of A as 
presented in Figurejsj Note, the coordinates (a;o, a:i, X2, Xa) are related to {x, u, y, v) 
by a linear change of coordinates. 

The above theorem implies that the set A — f]^yQP'^{B) is compact and con- 
nected. The following theorem implies that it is a non-trivial continuum. 

Theorem 2. The set A contains at least one fixed point and one period-two orbit 
for P. 
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In the literature there are already algorithms for computer assisted verification 
of uniform hyperbolicity and enclosure of attractors. The very pioneer and famous 
work has been done by Warwick Tucker jT . He proved that the Lorenz system 
for classical parameter values satisfies the geometric model by Guckenheimer and 
Holmes |GH) . Moreover, he gave the proof of the existence of the SRB measure 
supported on the attractor. 

Recently, Hruska |Hrll IHr2j proposed a method for computer assisted verifica- 
tion of uniform hyperbolicity. She successfully applied the method to the complex 
Henon map. The method proposed by Hruska uses the notion of box-hyperbolicity. 
This method reduces the verification of the hyperbolicity of an invariant set to the 
verification if some quadratic forms are positive definite. The box-hyperbolicity 
requires some conditions for the derivative of a map under consideration and its 
inverse. This is not a limitation in the case of the Henon map but makes the 
algorithm difficult to apply for ODE's. To compute the derivative of the inverse 
of a Poincare map we can either invert an interval matrix that is an enclosure of 
the derivatives of the Poincare map or integrate the ODE backwards together with 
the variational equations. In the first case the derivatives of Poincare maps are 
often non-invertible as interval matrices because of unavoidable over-estimations 
when integrating the variational equations. Backward integration does not help 
in many cases. After changing the time t — > — i in the equation the system often 
becomes stiff and again computationally difficult. In particular, this happens when 
the invariant set is an attractor with strong dissipation. 

Later, Aral [X] proposed a method for verification that a chain recurrent set 
for a map is uniformly hyperbolic. He applied the method to the Henon map [H] 
and verified that it is hyperbolic on the non-wandering set for a wide range of 
parameter values. This method, however, requires huge memory to work and it is 
computationally expensive by its very construction. Therefore, there is a little hope 
to apply it successfully in the high dimensional space and for nontrivial ODE's. 

Similar results were obtained by Mazur, Tabor and Kosciclniak jMTK| and by 
Mazur and Tabor fMT| . The authors introduced a notion of semi- hyperbolicity and 
applied the method to the real Henon map. 

Our method is similar in the spirit to that proposed by Hruska. We use a notion 
of strong hyperbolicity proposed in |KWZ1 [Z] as a theoretical tool for verifying 
uniform hyperbolicity. The main difference of the method proposed here is that it 
does not involve conditions for the inverse of the map under consideration. 

We believe that a successful application of our algorithms to a four dimensional 
Poincare map is a serious test for the method as well as the implementation and it 
proves their applicability. We would like to mention here that verification of uniform 
hyperbolicity for explicit given maps is significantly easier than for the continuous 
systems. It is clear that in a numerical approach to hyperbolic dynamics one needs 
to enclose both values and derivatives of the map under consideration. This is an 
easy task when the map is given explicitly, while for Poincare maps one needs to 
integrate the system with its variational equations. 

For rigorous integration of the system and its partial derivatives we used C° 
and solvers from the CAPD library |CAPD) that the author is one of the main 
developers. 

The paper is organized as follows. In Section [2] we give the notion of strong 
hyperbolicity and the theoretical background on how to verify uniform hyperbolicity 
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using this tool. In Section[3]we present algorithms used to define stable and unstable 
cones over the attractor. In Section[4]we give an application of the proposed method 
to the Smale map. We also present proofs of Theorem [T] and Theorem [2j 

2. Theoretical background. 

Let /: M" M" be a diffeomorphism and let AI C M" be a compact invariant 
set for /. We denote by TAI the restriction of the tangent bundle TM" to M. 

Definition 1. / is uniformly hyperbolic on AI if TAI splits into a direct sum 
TAI ~ © E'^ of two T f -invariant sub-bundles and there are constants c > and 
< A < 1 such that 

\\Df"\E^v\\<cX^v\\ and ||7?r"U.H| < cA"||«|| 

hold for n > 0. 

We will recall and reformulate the definition of strong hyperbolicity introduced 
in |KWZj and later extended in [Z]. 

Let = UiLi where Mi are compact sets having pairwise disjoint interiors. 
In our algorithms these sets will be boxes in some coordinate systems. Let us fix 
nonegative integers u, s such that u + s ~ n. Assume that at each set Ali we have 
fixed a linear coordinate system Ci. We define a quadratic form 

(4) Q{x,y) = \\xf-\\y\\^ xeR^,yeR\ 

For « = 1, . . . , iV we define positive and negative cones by 

Q+ ^{ue M" : Q{C,u) > 0} , 
Q- = {s e M" : QiCs) < 0} . 

We will denote this family by = (Q, {(Af^, C,;)},fl]^) and call cubical set with 
cones. 

Let A4 = {Q, {(Mi, Ci)}f^^ be a cubical set with cones and put AI ~ Uili 
For J = 1 . . . , we put = CjfC-\ 

Definition 2. / is strongly hyperbolic on M. — (Q, {(M^, Ci)}^]^) if for z £ Ali 
and J = 1, . . . , such that f{Mi) H Mj ^ the matrix 

[C,Df{z)CiYQ[CjDf{z)Ci'] - Q 

is positive definite. 

By Inv(/, AI) we will denote the maximal invariant set for / in A/, i.e. 
Inv(/, AI) ^{x£ AI : /"(x) S Af, /""(x) £ AI, for n > 0} . 

Theorem 3. // / is strongly hyperbolic on M. = (Q, {(Af^, Ci)}fl]^) then f is 
uniformly hyperbolic onT-L = Inv(/, M). 

The proof of the above theorem consists of several lemmas. 

Lemma 4. Assume f is strongly hyperbolic on M. = (Q, {(A/i, Ci)}f^]^) . 

(1) If f{Mi) n AIj ^ then for z £ AI, we have Df{z)Q+ C g+. 

(2) For z £ f{AIi) n Mj we have {Df{z))-^Q- C Q.^. 
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Proof. We will prove the first assertion. Since u £ we have Q[Ciu) > 0. Fix j 
such that f{Mi) n Mj 7^ 0. Put u = CiU and z = C^z. Since Q and Cj are linear 
and / is strongly hyperbolic on Ai we have 

Q{CjDf{z)u) = Q {DUj{z)u) > Q{u) = Q{C,u) > 0. 

Thus, by definition Df{z)u £ g+. 

We will now prove the second assertion. Fix z e f{Mi) D Mj and put z — 
f^^{z) G Mi. Assume that {Df{z))^^Qj <f. Q^, hence there is s e such that 
s :— {Df(z))^^s ^ . Reasoning as in the proof of the first assertion we conclude 
that 

Q{CjDf{z)s) > Q{C,s) > 0. 
Hence, s = Df{z)s G which contradicts the choice of s. □ 

Lemma 5. Assume f is strongly hyperbolic on M. = (Q, {(Afi, C^)}^]^) . There is 
e > such that for A e (f - e, 1 + e), i, j = f , . . . , such that f{Mi) n M, 7^ and 
z G Mi we have 

(5) [C,Df{z)C-YQ[C,Df{z)C-^] - XQ. 

Proof. The set of positive definite matrices is open. Thus for fixed i, j and z e Mi 
we can find a neighborhood Vij,z C M" x M of (z, 1) such that for (z, A) £ Vij^z the 
matrix 

[C,Dfiz)CiYQ[CjDfiz)Cr^] ~ XQ 
is positive definite. Then the assertion follows from the compactness of AI. □ 

Before we state the next lemmas let us define some constants. Put 

(6) Di_{M) = max{||CJ : i = 1,. .., A^}, 

(7) D2{M) = max{||Cr'|| : « = 1,.-.,A}. 

Lemma 6. Let A4 — (Q, {{Mi, Ci)}fLi) . There is a constant R > such that for 
u e M" and i = 1, . . . ,N we have \\u\\^ > R\Q{Ciu)\. 

Proof. We have 

\\uf = \\c-'c,uf>\\a\\-'\\cM\'> 

\\Q\\-^\\Tr,CM\^ > \m\-^Q{au) > D,{M)-^Qiau). 

Similarly 

> \\a\\-^\\7:yCM\^ > -D,iM)-^Qiau). 
Hence, the assertion follows with R — Di{Ai)^^. □ 

Lemma 7. Assume f is strongly hyperbolic on M. = (Q, {(Af^, Ci)}f^]^) . There 
are constants A > 1, c > such that for all i = 1, . . . , N , z € Mi D Inv(/, M), 
u e and k > we have \\Df''{z)u\\ > cA'^||u||. 

Proof. From Lemma [s] there is a constant A > 1 such that the matrices 

(8) y,,,- . = [C,Df{z)CrYQ[C,Df{z)Cr'] ~ XQ 

are positive definite for i, j = 1, . . . , A^ such that /(Afj) n Mj ^ and z e Mi. Since 
M is compact there is L > such that 

(9) t/,,,- ,H > ueM" 
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for i, j = 1, . . . , iV such that /(Mi) n Mj ^ and z e Af^. 

Fix z = Inv(/,M). Put Zk ^ f''{z) and Ak = Df{zk) for fc > 0. Let us fix a 
sequence {ik}k>o such that z^ G Mi^. Note, this sequence might be not unique. 
Put Ak = Ci,_^iAfcC~\ /c > 0. From ([8| for u 7^ and /fc > 1 we have 

(10) Q(C,,Df\z^)u) = Q[Ak-i ■ ■ ■ AoC.u) > 

\Q{Ak^2 ■ ■ ■ Af,a,u) >■■■> A'='iQ(ioC,„u). 

From Lemma [6] we get that there is a constant R > depending on Ai only, 
such that for m G M" holds 

(11) WDfizoM' > R\Q{C,,Df{zo)u)\. 
From (jSllol for u G (5+ we have 



(12) Q(ioC,„u) > Qia,u) + L\\a„uf > L\\a,uf > 

L\\C-'\\-^ur > LD^iM)- 



where D2{A4) was defined in (j7|. Combining (10 12) we obtain that for u G Q,^ 
and fc > 1 holds 

\\Df{zo)u\\>cX''\\u\\ 

with A = A^/^ and c — (RL)^/'^ {XD2{M))~^ . Clearly, adjusting the constant c we 
can obtain required inequality for fc > 0. □ 

There remains for us to show backward expansion in the negative cones. 

Lemma 8. Assume f is strongly hyperbolic on M. — (Q, {(M^, Ci)}f^]^) . There 
are constants A > 1, c > such that for z = 1, . . . , iV, z G Inv(/, M) n Mi, s G 
and k > we have ||£'/^'''(z)s|| > cA'^||s||. 

Proof. From Lemma [s] there is A G (0, 1) such that the matrices 

(13) y,,,- , = [CjDf{z)CrYQ[CjDf{z)Cr^] - ~XQ 

are positive definite for i,j — 1, . . . , N such that f{Mi)C] Mj ^ and z G Mj. Since 
M is compact there is L > such that 

(14) V,,jA^) > L\\4^ sgM" 
for i, j = 1, . . . , iV such that /(M,) n Mj ^ and z G M,. 

Let us fix z G Inv(/, M). Choose a sequence {ik}k>o such that Zk — f^''{z) G 
M,^ for fc > 0. For fc > put A^ = Df{zk), i/c = C,,_,AkC-^ and Sk = 
Df-^{zo)s. From ( 13 ) for fc > and s 7^ we obtain 

Q{C,,Df'''\zk)sk) = Q{A2 ■ ■ ■ AkC.su) > A'=-ig(Q,Sfc), 

that after substituting Sk and Zk becomes 

(15) - Q{C,,Df-''{zo)s) > -A-'=+iQ(C,,si). 

From Lemma [6] there is a constant R > depending on A4 only, such that for 
s G M" holds 

(16) ||i?r^(zo)sf > -RQ{a,Df-''{zo)s). 



From ( 13 14 ) we have 

Q(C,„s) = Q{C,,Dfiz,)si) = Q(iiC,,si) > Q(C,,si) + 
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Hence, for s e Q^^ holds 

(17) - Q(a,si) > L\\C\,4^ - Q{a,s) > L\\C\,4' > LD2{M)-^\\sf. 



Combining (15 17) we obtain that for fc > 1 and s e Q^^^ holds 

\\Dr\zo)s\\>c\''\\s\\, 

with A = > 1 and c = (Li?)^/^ (XD2{M)y^ . Adjusting the constant c if 

necessary we can obtain the same inequality for fc > 0. □ 

Proof of Theorem^ Let us associate to each z € Inv(/, M) an index iz such that 
z £ M-i^ . This index might be not unique but for each z we can fix one of them. 
For simplicity we will write M^, C^, Qf instead of M^^ , Ci^ , Qf . 

From Lemma [4j Lemma [7] and Lemma |8] there are constants c > and A > 1 
such that for z £ Inv(/, AI) we have 

Df{z)Qt c g+ 



Dfiz)-'Q- C Q-,, 



\\Df''iz)u\\ > cX''\\u\\, forueQ+, 
\\Df-''{z)s\\ > cA'^PII, forseQ;. 

Note, the cones Q+ and Q~ are disjoint and the sum of them spans M". Hence, 
the assertion follows from |KH| Corollary 6.4.8]. □ 



3. Algorithms. 

3.1. Graph representations of maps. Let X be a topological space and let 
N C X. For a family of subsets U E 2^ by supp (U) we will denote its geometrical 
representation supp {U) = iJueu ^■ 

Definition 3. We will say that lA E 2^ is a cover of N if lA is a finite set and 
N C supp(W). 

Definition 4. LetlA,V he a finite sets. A multivalued function f : U V is called 
a combinatorial map. 

Covers of compact sets as well as combinatorial maps appear naturally in the 
computer assisted proofs in dynamical systems. We usually cannot compute the 
range of a map f : X Y over a domain U C X. Instead, one can fix covers X, y 
of X and Y, respectively, and then compute combinatorial representation of the 
map in the covers X,y. This is usually realized in the following steps. 

• For U d X compute a minimal cover 

U ^f]{V (Z X -.U C supp (V)} 

• For Vi €U compute a cover Vi C 3^ of f{Vi). This step is often realized by 
means of the interval arithmetic [M]- In most cases it is very difficult to 
compute Vi as a minimal cover of J{Vi). 

• Define a cover of f{U) by IJ V^. 

The above considerations lead to a very natural definition of the combinatorial 
representation of a map. 
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Definition 5. Let f: X Y be a map and let us fix covers X,y of X and Y, 

respectively. 

We will say that a combinatorial map J- : X —o y is a combinatorial representa- 
tion of f if for X € X and for allU ^ X such that x € U holds f{x) G supp (J-"(C/)). 

For a map f : X ^ X and a cover A" of X it is natural to encode a combinatorial 
representation of / as a directed graph. A finite directed graph is a pair Q = (V,£), 
where V is a finite set whose elements are called vertexes, and £ C V x V is a set 
of selected (ordered) pairs of vertexes. The elements of £ are called edges. 

Definition 6. We say that a graph Q ~ (V, £) is a representation of a combinatorial 
map J-:U^U ifV=lA and 

£ = {([/, V) (^U xU -.V ^{U)}. 

We will need a notion of the transposed graph and outgoing edges. For a graph 
Q = iy,£) by — (V',£') we will denote a graph in which V' = V and 

£' ^{{U,W):{W,U)(.£}. 

For a vertex F in a graph Q = (V,f) by out{V,Q) we will denote the set of 
vertexes 

out{V,g) = {P eV : {V,P) e £}. 

3.2. Enclosure of an attractor. In this section we will describe an algorithm that 
we used to enclose an attractor. There are already existing packages for rigorous 
enclosures of invariant objects, see for example |GAIOj . We used, however, our 
own implementation dedicated for this purpose which seems to be easier and makes 
the software independent on the other libraries. We will use the fact that what we 
want to enclose is an attractor. 

The idea is very simple. First, we fix a cover X of the observed attracting 
domain. Then we choose a box V from observed attracting region, and we enclose 
its forward trajectory using the sets from the cover X as long as the trajectory 
does not leave supp {X). Since the cover X is finite by its definition, after a finite 
number of steps there are no new sets in the cover of the trajectory of V and we 
can stop the procedure. This approach generates an enclosure of some invariant 
set, that is expected to be an attractor. 

The Algorithm [T] summarizes the above considerations. 

Lemma 9. Fix a cover X of X. Assume Algorithm^ is called with its arguments 
V , Q — (V,£) and f . If the algorithms stops without the Failure result then 

(i) the combinatorial map J-": V — o V encoded by the graph Q — (V,£) is com- 
binatorial representation o//|supp(v)) with the natural cover V o/supp(V), 

(ii) supp (V) is a positive invariant set for f, i.e. /(supp (V)) C supp (V), 

(iii) U::or(^)csupp(V). 

Proof. First observe that the algorithm always stops. By contradiction, assume it 
does not stop. This means that in each iteration of the main repeat-until loop (lines 
4-15) the set U is not empty. This implies, that the number of elements in the set 
V, which is updated in line 14, is strictly increasing. But this is a contradiction 
with V C X, which is a finite set by the definition of cover. 

Let G be the result returned by the algorithm and let : V ^ V be the combi- 
natorial map represented by the graph Q. 
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Algorithm 1: enclose forward trajectory 



Input: V : element of X, 
g^iV,£) : graph, 
/ : a map 

Result: Q - graph 

Data: U,J- : subsets of X, 
G' = (V',£') : graph 

1 begin 



2 


(V,£ 


)^(0,0); 


3 


U ^ 




4 


repeat 


5 






6 




foreach U eU do 


7 






if f{U) (f. supp (A") then 


8 






1^ return Failure; 


9 






T ^ [/([/)] at; // computed cover of /(J7) 


10 






V ^ V U {[/} U J"; 


11 






foreach F <E T do 


12 






L £:'^£:'u{(c/,F)}; 


13 




U \ {VUU); 


14 




(V,£)^(V,£)U(V',r); 


15 


until = ; 


16 


return Q; 


17 end 







We are proving the assertion (i). Let x G supp (V) and let e V be such 
that X G V. Observe, that if J-{V) is not empty then f{V) C supp {J- {V)) - as 
guaranteed in lines 9, 11, 12 and 14 of the algorithm. Hence, in this case the 
assertion follows. To finish this step we observe, that the main repeat-until loop 
stops when the set U is empty. But this set contains of exactly these sets U E V 
(lines 10 and 13) for which the value J^{U) is not computed yet. Hence, after the 
algorithm stops, for each V GV, supp {out{V, G)) contains a computed enclosure of 

fiv). 

Assertions (ii) follows from (i) and the definition of combinatorial representation. 
To prove (iii) it is enough to observe that V C supp (V) as guaranteed in line 3, 
10 and 14 of the algorithm. Then the assertion follows from (ii). □ 

Using Algorithm [l] one can compute combinatorial representation of a map re- 
stricted to some positive invariant set, which is expected to be an attractor. For 
further applications this graph should be as small as possible and it should still 
contain enclosure of the invariant set of / restricted to supp(V). Therefore it is 
important to choose the starting box V such that it contains points from the attrac- 
tor. An easy way to do that is to take a point from the observed basin of attraction 
and compute non-rigorously its forward trajectory for some long time. Then we 
can choose a set V from the cover that contains the last computed point. 
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3.3. Computing of coordinate systems. In this section we assume that the 
graph Q = (V, £) is a combinatorial representation of / resulting from Algorithm[l] 
We present a heuristic method for computing of coordinate systems Cy at the 
vertex V, ior V G V. We assume that at the beginning the coordinate systems are 
not set, which we will denote by Cy=NULL. 
The method consists of the following steps. 

Step 1. Finding cycles in the graph. First we try to find (non-rigorously) as 
much as possible of periodic points of / in supp(V). This is important, because 
given a period p point, say x, a very natural choice of a coordinate system in a vertex 
containing this point is the Jordan basis of DfP{x). Of course it may happen that 
several periodic points belong to the same vertex. In this case we are choosing one 
from these points of the lowest principal period. 

More formally; let us fix a positive integer N. This is an upper bound for the 
highest period of orbits we will search for. Since Q is a combinatorial representation 
of /, all periodic points of / must belong to some cycles in the graph Q. Let Vi, 
i = 1, . . . N he a. set of vertexes V such that there exists an i-periodic cycle in the 
graph Q through V and no lower period cycle containing V exists. The sets Vi can 
be computed by means of standard graph algorithms and we omit details here. 

Step 2. Refining cycles to periodic points. In this step we compute the sets 
of points Pi, i = \,...,N that are good approximations of periodic points. For 

V ^Vi we proceed as follows. Let x be an approximate center of V . The point x is 
a seed point for the standard Newton method for zero finding problem applied to 
the map g{x) = p{x) — x. If the Newton method converges to a point y e supp (V) 
which is of principal period i then we insert the point y to the set Pi. 

We would like to comment here, that we do not try to find all periodic orbits of 
low periods for the map /. Using our approach one can detect most of them if the 
set V consists of small enough sets. Clearly, a larger number of detected periodic 
points is better for setting the coordinate systems but the task of finding all low 
period orbits is not a goal for us. 

Step 3. Setting coordinate systems in a neighborhood of periodic points. 

The goal of this step is to assign a coordinate system to vertexes that contain 
computed approximate periodic point. For y d Pi we proceed as follows. Let 

V G V he a vertex such that y G V. Such a vertex exists since we know that 
y G supp (V) - as guaranteed by the previous step. 

• If the vertex V has already assigned a coordinate system then we proceed 
with the next point from Pi. 

• We compute an approximate derivative Df^{y) and an approximate matrix 
My of normalized eigenvectors. The columns of My are sorted by decreasing 
absolute value of the associated eigenvalue. 

• We assign to the vertex V a coordinate system Vc = My. 

We would like to comment here that the sets Pi are proceeded by increasing index 
i. This guarantees that a vertex V will have assigned coordinate system from some 
periodic point of the lowest principal period. Note, the resulting coordinates depend 
on the order chosen in the set Pi. 
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Step 4. Spreading coordinates from periodic points. The main idea is 
to propagate the coordinate systems from periodic points by the action of the 
derivative. It is well known that the direct propagation of coordinates by Df 
in floating point arithmetic usually results in collapsing of these coordinates to a 
singular matrix. Therefore we will involve some orthonormalization process. 
The Algorithm [2] describes the procedure. 



Algorithm 2: spread coordinates from periodic points 
Input: g = {V,£) : graph, 
<S : subset of V, 
fc > 1 : integer 

Result: {Cv}vev - coordinate systems at each vertex in the graph Q 
Data: u : vector; 

M : stack of float matrices; 

C : matrix; 

V : subset of V; 

i : integer 



begin 
repeat 

V ^ 



foreach F G 5 do 

u centre(y); 

/ / we will save on the stack derivatives along the trajectory; 
C^Cv; 
for i — 1 to k do 
push(i^/(u),M); 
C ^ Df{u) ■ C; 
u ^ f{u); 

1 1 then we compute coordinate system at the image; 
C -s— GrammSchmidtOrthonormalization(C); 
for i — \ to k — 1 do 

(top(A/))"^ • C; 
pop(M); 

pop(Af); //to free the memory; 
C -s— normalizeColumns(C); 
foreach W G out(y, G) do 
if Cw ^NULL then 

Cw C*"^; 



S 

until S = ( 



25 end 
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Lemma 10. Assume the Algorithm^is called with its arguments Q = (V,f ), 5 7^ 
and k > 1. Assume Cy 7^ NULL for V (Iz S. If Q consists of only one strongly 
connected component then the algorithm stops and Cy 7^ NULL for allV £V. 

Proof. We will prove that the algorithm always stops. Assume it is not the case. 
Then in each iteration of the main loop the set V is non-empty. But this set contains 
(lines 3 and 19-22) exactly these vertexes for which the coordinate system is assigned 
in the current iteration of the loop. Hence, the set of vertexes at which we have 
assigned the coordinate system will increase in each iteration. This contradicts the 
fact, that V is a finite set. 

We will prove that the algorithm sets coordinates in each vertex. Assume it is 
not the case. Take Q"^ = (V',£')- Let be a vertex such that Cv=NULL. Then 
for U e out(l/, V) holds Cu = NULL. Otherwise, U wih belong to the set S in 
some iteration and then the coordinates will be set in V . Denote Vi = out(V, G'^)- 

In a recurrent way we can define 

Vfc-|J{out(C/,g^) :C/e Vfe_i}. 

The same reasoning proves that Cc/=NULL for U ^ Vk- 

Since G consists of only one strongly connected component, G^ does too. Hence, 
for some if > we have Vk = V which contradicts the assumption that Cu ^^NULL 
for ;7 G 5 and 5 7^ 0. □ 

The input argument S is the set of vertexes that contain periodic points and for 
which the coordinate systems were already computed in Step 3. As we have seen, 
lines 5-18 are not important for the correctness of the Algorithm[2j In these lines a 
heuristic method for computing of coordinate systems is proposed. We propagate 
the actual coordinate system Cy by the derivative along a short part (fc iterations) 
of the trajectory of the center of V. Then we perform orthonormalization of columns 
and propagate the obtained matrix backwards k — 1 times. This heuristic can work 
because for the inverse map the stable directions are attracting for Tf and after 
few steps we get some approximation of stable directions at /(centre(y)). 

In this step we also assumed that the graph G consists of only one strongly 
connected component. It is not a restriction, since the method is dedicated for 
attractors. In practice, Algorithm [T] always returns such a graph if the initial set 
V is chosen carefully. Formally, it is easy to verify that a graph consists of one 
strongly connected component using for example Tarjan's algorithm |Ta) . 

Remark 11. The result of the algorithms described in steps 2,3,4 depend on the 
order chosen to index the sets of vertexes and points (Pi, Pi and S). Moreover, 
these algorithms are not deterministic, when implemented in parallel mode. The 
coordinate system at each vertex is computed by the first thread that reaches this 
vertex. 

Anyway, the conclusion of Lemma \10\ holds; we have computed some coordinate 
.system at each vertex and we can start verification of the cone conditions. 

3.4. Verification of the cone conditions. This step is quite straightforward. 
Let G = ) be a graph which encodes a combinatorial representation of a map 
/. Assume we have computed some coordinate systems {Cv}vgv- The goal is to 
check whether the assumptions of Theorem |3] hold true for /. This can be done by 
simply verifying the positive definiteness of some interval matrices corresponding 
to every edge in £. 
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Algorithm 3: verify cone conditions 

Input: g = (V,f) : graph; 

Q : matrix; 

/ : map; 
Result: 14 : subset of V; 
Data: V, W : element of V; 

D, M, A : interval matrix; 



begin 



foreach V E V do 

D ^ [Df{V)]; II computed interval enclosure of Df(V); 
foreach W e out{V,g) do 
M ^ [CwDCy^]; 

-Q-M-Q; 
if cannot verify that A is positive definite then 
[_ U=UU{V}; 



return U; 



11 end 



We have the following obvious lemma. 

Lemma 12. Assume the Algorithm^is called with its arguments Q — {V,£), Q 
and f . Assume that all the matrices {Cv}vi£V olt^ invertible. The algorithm always 
stops and returns a set U such that for V Cz V \ U and W G out(V, G) the interval 
matrix 

[CwDf{V)CyYQ[CwDfiV)Cy'] - Q 

is positive definite. 

A direct consequence of the above Lemma and Theorem [3] is the following 



Corollary 13. The same assumption as in Lemma 12 Assume Q encodes a com- 
binatorial representation of the map f and let Q be quadratic form as in If the 
Algorithm^ called with its arguments Q — (V,£), Q and f returns an empty setlA 
then f is uniformly hyperbolic on Inv(/, supp (V)). 

4. Applications. 

In this section we will show how the method introduced in the last section works 
in a very easy example. Then wc will give a proof of Theorem [ij 

4.1. Toy example - Smale map. A natural choice to test the algorithms and 
implementation is a very simple example. Here we have chosen the well known 
Smale map s : x [0, 1] ^ ^ [o, 1] 

(18) s{x,y,t) = (0.1a; + 0.5cos(27rt),0.1?/ + 0.5sin(27rt),2i mod 1). 

The map is proved to be uniformly hyperbolic just by simple hand-made calcula- 
tions. It is easy to see that the set [—1,1] x [—1,1] x [0, 1] is positive invariant for 
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s. Define a cover of this attracting domain of the form 

We will call the number k in Xf; the resolution. After application of Algorithm [T] we 
got an enclosure of the attractor consisting of 567 boxes for fc = 4 and 3333 boxes 
for k = 6. These enclosures are shown in Fig. [2] 




Figure 2. Enclosure of the Smale attractor for the resolutions 
k — A and k — Note, the planes t = and t — \ should be 
identified. 



We completed the proof of uniform hyperbolicity for the resolution k = A. We 
run the algorithm for finding cycles in the graph and refining them to periodic 
points with the maximum period set to iV = 3. They returned 1 candidate for fixed 
point, a period 2 orbit and two period 3 orbits. 

Then we applied the algorithm for spreading coordinate systems with k = 2 
(this is the number of forward iterates when propagate the coordinate systems) 
and verification of the cone conditions. We run Algorithm |3] with the quadratic 
form 





-1 







-1 



and it returned an empty set lA of unverified vertexes. The program executes within 
less than 1 second on a laptop-type computer. 



4.2. Application to the system of coupled van der Pol oscillators proof 
of Theorem [ij Consider the system ([!]) with the parameter values (|3| . Let P be 
the Poincare map as defined in ([2]). Let us perform a linear change of coordinates 



Xq xjv^^ X\ (w ^ux)l'^li '^2 y C.yxX Cyy^U^ X3 V Cyx"^ ^vu ^vyV: 
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algorithm 


wall time (h:mm) 


comments 


enclosure of attractor, 
Algorithm [l] 


2:16 on 224 CPUs 


7 970 392 boxes 


finding of cycles in graph, 
max period 6 


0:58 on 32 CPUs 


2190 cycles found 


finding of periodic points, 
max period 6 


0:06 on 32 CPUs 


105 points found 


computing of coordinate systems. 
Step 4 and Algorithm [2] 


6:59 on 32 CPUs 


parameter k = 2 


verification of cone condition. 
Algorithm Is] 


4:24 on 224 CPUs 


empty set 
of unverified boxes 



Table 1. Setting, results and time of computation of the algo- 
rithms applied to the map P. 



where 



= 0.438, 
Cyu = 0.226, 

= 0.029, 
ro = 0.812, 



-0.042, 



= -0.218, 



Cy 

C,y = -0.118, 

ri = 0.721. 



As it is observed in |KSe| . in these coordinates the attractor is more aligned to the 
coordinate axes. In the sequel we will consider the Poincare map P = L^^PL, 
where L is the matrix of coordinate change. 





Figure 3. Projection of the enclosure of the attractor for P com- 
puted with the resolutions k = 14. 



After few attempts we found that the resolution fc = 14 is enough for verification 
of the hyperbolicity of the observed attractor for P. This means that we have used 
a uniform grid in some bounded domain of the form 

A- . I [-^- + ^\-[b,b + l]>^[c.c+l]^[d.d+l] ^^^^^^^^^ 2,^, _ ^ 

In Table [T] we give the results of applying of the algorithms for generating enclo- 
sure, setting of coordinate systems and verification of strong hyperbolicity. We run 
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Algorithm [3] with the quadratic form 






-1 









-1 



and it returned an empty set of unverified vertexes. Hence, the uniform hyperbol- 
icity of the set Inv(-P, supp (V)) is verified. To complete the proof of Theorem [l] 
we used the CHOMP library |CHOMP] in order to compute homology groups of 
the cubical set represented by V. The program verified that Hq{V) = Z, hence 
the set supp (V) is connected. The group Hi{V) also was nonzero, hence the set is 
non-contractible. 

This finishes the proof of Theorem [T] 

4.3. The existence of a fixed point and a period two orbit proof of 
Theorem [2j As a byproduct of the proof of Theorem [T] we got a good numer- 
ical approximations of over one hundred periodic points. To prove that the set 
Inv(P, supp (V)) is nontrivial we used the interval Newton operator [M] as a tool 
for verification of the existence and local uniqueness of zeros of maps. We have the 
following lemma. 

Lemma 14. Assume f is a smooth map in a neighborhood of an interval set X C 
M" and fix x G intX . If the interval Newton operator 

N{x, XJ)^x- [Df{X)]Y\X ~x)c mt{X) 

then the map f has unique zero in the set X. Moreover, if x^, d X is the unique 
zero of f in X then a;* € N{x, X, f). 

Proof of Theorem [1[ A good numerical approximation of a fixed point for P is 

-0.2217320903523841312 ^ 
0.97492087818262221051 
0.018539814418960732759 
-0.03127925657964333925 

— this is a point returned from the step of finding periodic points in supp (V). Let 
X be a ball of radius lO^^'' in the maximum norm centered at x. We computed 
the interval Newton operator N{x, X,P — Id) and we got that 

\\N{x,X,P -Id) -5||o < 2.5- 10"^^ 

Hence, N{x,X,P — Id) C mt{X) and the map P — Id has a unique zero in X. 
Clearly this is a fixed point for P. We also verified that N{x, X,P — Id) C supp (V) 
and this finishes the proof that P has a fixed point supp (V). 

In a similar way we verified that the attractor contains a period two orbit. It is 
enough to show that the map F{x,y) = {P{x) — y,P{y) — x) has a zero for some 
{x.^,,y.^,) and x.^, y*. We have a good candidate for periodic points 



1.0708168699357167158 
0.093670228872794565237 
0.015571609254819055212 
0.024616361886390349154 



-0.22099700889399603573 
-0.706209321733317795 
-0.044189820507618084004 
-0.006451355101601539937 
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Let Z C be a ball of radius 10"^° in the maximum norm centered at z = (x,y). 
We verified that 

\\N{z, Z,F)-z|lo<2.10-i2. 

This proves that F has a unique zero G Z. Since ||a; — j/|lo > 2 • 10"-^° 

we have ^ and both x^,y^ are period two points for P. We also verified 
that {x^,y^,) e N{z,Z,F) C supp (V) x supp (V) and this completes the proof of 
Theorem El □ 



5. Conclusions. 

In this paper we proposed a method for rigorous verification of uniform hyper- 
bolicity for maps. The method has been successfully applied to a Poincare map of 
a non-autonomous ODE in four dimensions. We believe the method can be applied 
to other systems. 

Kuznetsov [K2j proposed also a non-autonomous system on the plane which ap- 
parently possesses an attractor of Plykin type in the Poincare map. The equations 
are quite complicated and given by 
(19) 

i — — 2£y^rii(a;, y, <) (cos cos — a; sin (l^ cos 1^)) -|- 
Ky^2{x,y,t) (cos (^j sin |i) — a; sin (^ sin sin ^t, 

y = £j/fii(x,?;,i) (2a;cos (f cosfi) + (1 - x2 + y2)sin(|cosfi)) - 

Kn2{x,y,t) (a;cos(f sinft) + i(l - + ?;2) gin (f sin f t)) sin f i, 

where 

2x008 (f cos ft) + (1 - a;2 _ y'^) gin (| cos ft) 



ni{x,y,t) = 
^2{x,y,t) = 



(1 + a;2 + j/2 

-2a; sin (f sin ft) + (1 - a;^ _ y2) cos (f sin f t) 1 

^2 ' 



(l + a;2+y2)2 ^2' 

Kuznetsov jK2] did a numerical study of this system and observed that for a wide 
range of parameter values around K = 1.9 and e = 0.72 the Poincare map has a 
hyperbolic attractor. A rigorous enclosure of the attractor for K = 1.9, e = 0.72 
resulting from Algorithm [l] with the resolution A: = 9 is shown in Fig. [4j The 
resolution k = 9 was not enough to perform the step of verification of uniform 
hyperbolicity of this map. We also enclosed this attractor with resolution fc = 12 
but again this was not enough to verify the hyperbolicity. Further subdivision was 
not possible due to memory limitations in the computer that we used. 

Computations for the system (19 1 are much harder than for ([I]). First - the 



vector field is a quite large expression and even if the system is lower dimensional, 
it is more complicated for rigorous integration. Second, and most important; as 
easily seen in Fig.]!] the attractor is very thick and harder to cover than the almost 
one dimensional Smale solenoid. An approximate Hausdorff dimension reported in 
jK2j is ca 1.7. We reached the limit of available memory on our computer when 
enclosing the attractor with a thinner cover. We believe, however, that some further 
optimizations in data structures will help to prove hyperbolicity of this system. A 
straightforward idea is to use non-uniform covers. 
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Figure 4. Enclosure of the attractor of the Poincare map for the 
system ( 19 ) with parameter values K — 1.9, e = 0.72 and with the 
resolution fc = 9. 



5.1. General invariant sets. In this paper we verified hyperbolicity of attractors. 
But the method can be easily extended to invariant sets. The only part which 
requires some modification is the method for generating an enclosure of an invariant 
set. If the invariant set is an attractor, the strategy of "inner" enclosure is better 

- Algorithm [l] For a general invariant set an "outer" approximation is necessary. 
We also implemented "outer" approximation for enclosing of general invariant sets 
but we do not report details here, since it is rather standard procedure. As a test 
case we proved that the well known Henon map [Hj 

'Ha,b{x, y) = {l + y- ax"^, bx) 

is uniformly hyperbolic on the invariant set for the parameter values a — 5.4 and 
b^-l. 

As mentioned above, we used an "outer" strategy to enclose the invariant set 

- the result is shown in Fig. [5] and it consists of 8832 boxes. Then we applied 
algorithms for finding cycles and periodic points, computing of coordinate systems 
and verification of the cone condition with the quadratic form 



Q 





--1 



The program executes within less than 1 second on a laptop-type computer with 
the Intel Core 2 Duo 2GIIz processor. A computer assisted proof for the same 
parameter values is presented in [3 IMT| . In |MTj the authors report the time of 
computations less than 10 seconds on a comparable CPU. 



5.2. Implementation. The C++ program for verifying the hyperbolicity has been 
implemented by the author and the source code is available at |W]. The code is 
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Figure 5. An enclosure of the invariant set for 'Kaf, map with the 
parameter values a = 5.4, h = —1. 



highly parallelized using the Open MP library supported by the compilers gcc-4.2 
or newer. The code is written in a very generic way; it heavily uses template tech- 
niques. In fact, an input to the algorithms is an abstract map (template parameter) 
for which we assume that we know how to compute values and derivatives. In the 
case of Poincare maps we used for this purpose the and solvers from the 
[CAPDj library - see also |Z1| for an efhcient algorithm for integration of first order 
variational equations. 

All computations regarding Theorem[l]and Theorem [2] were performed on a clus- 
ter with 7 computers, each having 8 Quad-Core AMD Opteron(tm) 8354, 2.2GIIz 
processors with 64GB of RAM. During the computations the maximal usage of the 
memory was 32%. 
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